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VJ ABSTRACT 

Two problems are presented which are linear on two adjacent inter- 
vals but not on their union. These problems are associated with the 

(AX + BF, 0 <t <ti 

differential equation X = ^ t^<t<T > where X is the matrix 

[y) 9 where F is a 2 x 1 matrix, and where A,B,C, and D are 2x2 matrices 
of functions of t. t^ is a variable, hence the differential equation is 
non-linear. Problems associated with this differential equation are 
called semi-linear . 

In the first problem, a condition is found on t^ and F which must 
be satisfied whenever x(T) is to be a maximum with y(T) fixed. In the 
second problem, conditions on F and t^ are found which must be satisfied 
for x(T) to be a maximum for a fixed y(T) and for a fixed x(tx) • A numer- 
ical routine is developed which yields successive approximations to the 
maximum value of x(T). 

The basic theory of the methods is presented, and the problems are 



developed in the context of optimum control. 
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I. Introduction. 



Nonlinearities make the solution of optimum control problems more 
difficult. In general, different techniques must be developed for each 
new nonlinear problem. 

The particular problems to be discussed are those in which a set of 

functions x. of a variable t is related to a second set of functions f. 

1 1 

of t through a differential equation which is linear on each of the in- 
tervals (0,t^) and (t^,T). If t^ is a variable, the differential equa- 
tion, though linear on each of the intervals and (t^,T), is 

non-linear on the interval (0,T). Such a differential equation will be 
called a semi-linear equation; problems in which a semi-linear differen- 
tial equation occurs will be called semi-linear prdblems. Such problems 
arise in the theory of optimum control; the problems to be discussed 
will be in the terminology of optimum control theory- 

The specific semi-linear problems to be considered arise from the 

differential equation X = ^ 9 where X,A,B,C,D, and F 

CCX + DF , 1 2 < t < T 

are described as follows: X is the matrix (^) whose elements are con- 



tinuous functions of t; X is the matrix of derivatives of x and y 
with respect to t. A,B,C, and D are 2x2 matrices whose elements are 
functions of t which are piecewise continuous and bounded on the inter- 
val (0,T). In addition it is necessary that B and D be nonsingular. F 



0 O' 



is the matrix [ of functions of t; associated with F there is some 
■2 

constraint, specified in each problem. F is called the matrix of con- 
trol variables , and those matrices F meeting the specified constraints 
are called allowable . 

Associated with the above differential equation are certain boundary 
conditions which X(t) is to satisfy. Generally the initial point is 



5 



given. In addition, one or more components of X may be specified at some 
value of t on the interval (0,T). Curves satisfying the given boundary 
conditions and for which F is allowable, are called admissible . 

In the terms that have been defined above, and for the above differ- 
ential equation, the semi-linear problems to be solved are the following: 
(1) Find conditions on F and on t^ which must be satisfied if x(T) is to 
be a maximum for a given y(T) , where t^ must be determined and where x(t^) 
and y(t^) are unspecified. (2) Add the constraint that x(t^) is to be 
specified, for a value of t^ to be determined. Then (a) find conditions 
on F and t^ which must be satisfied for a curve to be admissible; (b) find 
additional , conditions on F and on t^ such that an admissible curve will 
yield the maximum x(T) ; and, (c) develop a method of successive approxima- 
tions which will give values of x(T) converging to this maximum. 

The following results are obtained: In the first problem, a neces- 

sary condition at tj^ is derived. In the second problem, conditions for 
a maximum are derived, and a numerical method is given for obtaining this 
maximum. These results are another step in the solution of nonlinear prob- 
lems of optimvim control. 

I wish to thank Professor Fpulkner for the encouragement and help he 
has given me and for the permission to use class notes of his course, 
Differential Equations of Optimum Control. 
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II. General terminology and some sufficient conditions for the linear 
problem. 

In this section we define the linear problem and explain the termi- 
nology to be used. We then derive some sufficient conditions for the 
linear problem. 

Let us consider the differential equation 
(1) AX + BF 



where X is the n x 1 matrix (x^) , where A and B are the n x n matrices 
(a^j) and , respectively, and where F is the n x 1 matrix (f ^) . We 

will always assume that B is nonsingular, although the results apply to 
many problems where B is singular. The the , and the f^ are 

functions of t that are bounded and piecewise continuous on the inter- 
val (0, T) . The a^j and the b^^ are givemif unctions ; the f^ are func- 
tions satisfying a given constraint but otherwise unspecified. The f^ 
are often called control variables and the x^ state , or dependent 
variables . 

Associated with the differential equation (1) are certain boundary 
conditions which we want X(t) to satisfy. Generally the initial point 
is given. In addition we may specify that certain elements of X take on 
stated values at various fixed points of the interval (0, T) . Note that 
the differential equation is linear , since A and B are functions of t 
only. 



If we multiply both sides of equation (1) on the left by the 1 x n 

T . i 

matrix K (the transpose of the matrix K) , whose elements k are func- 
tions which have not yet been specified, and integrate from t=0 to t=T, 
we get 



( 2 ) 



T 

^ k'^X dt 

0 



T 

^k’^AX dt 

0 



+ 



T 

^k'^BF dt. 
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Integrating the left side by parts and rearranging terms, we get 
T T T 

= + K^A)X dt + ^k’^BF dt. 

0 0 0 

Now let us choose and take the transpose; we get 



(3) 



T 

K^X 



(4) K = -A K. 



This is the adjoint equation associated with equation (1). Choosing K 
as a solution to the adjoint equation and substituting into equation (3) 
we get 

T T 



(5) 



T 

K X 






k'^BF dt. 



0 0 






we can choose a solution K to the adjoint equation so that K (T) is 
the matrix (10 0 ... 0). If we substitute this solution into equa- 
tion (5) and rearrange terms, we get a solution for x^(T), namely 

T 

(6) x^(T) = K^^(0)X(0) + ^ K^^(t)BF dt. 

0 

In the same way we chose K^, we can choose K , K , ... , K such that 
iT 

K (T) = ^^i^2i ^ni^ ’ where 6^^ equals zero, if i/j, or one, if 
i=j. Suppose now that we have the n linearly independent column mat- 
rices chosen above. Such a set of n linearly independent matrices 
is called a fundamental set of solutions for equation (4). If we form 
the n X n matrix 7^ whose i'th column is the matrix K^, then we may 
write every solution to equation (4) in the form where C is an 

n X 1 matrix of constants. Conversely, every product of the form/^C, 
being a linear combination of solutions to equation (4), is also a 
solution. Furthermore K is itself a solution. This is shown as fol- 
lows: Take any product }(^C, where C is an n x 1 matrix of constants 
and where ^ is the matrix defined above. This product is a solution 
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to equation (4). Hence = -A^T^C. But this equation is valid for 

every G and hence for the C whose transpose is (1 0 0 ...0). Substitu- 

ting this C into the above equation, we see that the first column of rC 

T 

is the same as the first column of -A But we could have equally well 

chosen the i’th element of C as one with the others zero; we hence could 
have found that the i*th columns of both sides were equal, for i=l, 2, 

s 

n. Hence is a solution to the adjoint equation, which 

is what we wanted to show. Furthermore, if 6^ is a small variation in 
^ , and if 6 is the corresponding variation in then /C ^ ^ ~ 

• ^ rn 

-A^( i-A- -a'^^/T. since K - Hence the varia- 

tions in K are related to the variations in 7C by the adjoint equation. 

We have shown that solutions to the adjoint equation can be used to 
find solutions to the equation X = AX + BF when F is known. The next 
problem we are concerned with is that of finding F so that for a given 
X(0) and a given T, a specified component of X(T) is a maximum. F can 
be regarded as a column matrix of forcing functions; these forcing func- 
tions are the components of the forcing function vector . There may be 
one of several types of constraints on F. In rocket thrust scheduling, 
for example, the acceleration possible at any one time is limited by a 
function of the mass of fuel on board at that time. If we call this 
function cp(t) , the corresponding constraint on F is that |F|6V(t) , 
where |f| is the square root of the sum of the squares of the elements 
of F. Another type of constraint on F is |f^|^|a^|, where a^ is a 
given function of t. Problems of this second kind are called bang- 
bang control problems^, % In every problem it is necessary to state 
the constraint on F; forcing functions which satisfy the stated con- 
straints will be called allowable. Solution curves to equation (1) 
for which F is allowable will be called allowable curves: allowable 
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curves satisfying the differential equation and satisfying the given 
boundary conditions will be called admissible * 

We want a method for choosing an allowable F such that X(0) is a 
given point and such that a specified element of X(T) is a maximum* A 
most important principle enables us to state conditions which F must 
satisfy whenever the desired maximization takes place, namely Po ntry- 
agin’s Maximum Principle ; The control variables must be chosen from the 
set of allowable controls so as to maximize a scalar product of some 
solution to the adjoint equation and the forcing function vector at every 
time t*:[x3^ will use this principle extensively in the following 
pages . 

Consider the following problem and see how the maximum principle 
may apply to it* Suppose that we want a curve beginning at some fixed 
point at time t = 0 on which some element of X at t=T, say x^(T), is a 
maximum and such that on the curve these three conditions hold: First, 

F is allowable* Second, the curve satisfies the differential equation 
X = AX + BF for all values of t between t=0 and t=T where T> 0 is given* 

Third, x^(T) = for i=2, *** , N, where x^^ are given, and where 

KN^i^n, and where Xj^^^(T), *** , are unspecified* N=1 means that 

no values of the x^ are given at t=T; N=n means that all values except 

x^ are given at t=T. A proof given by Faulkner l3l proves that the fol- 

lowing hypotheses are sufficient for a curve C* to yield a maximum x^(T). 
Suppose that we have found a curve C*^ with forcing function F and have 
found at the same time a solution K* to the adjoint equation which to- 
gether satisfy the following hypotheses: 

Hi. C* is admissible: it begins at X^ and ends on the manifold de- 
fined by X 2 (T) = X 2 .J,, ... , “ ^NT’ allowable* 
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H2. F* maximizes K*^BF for all allowable F, i.e. K*^BF*^ K*^BF for 
all allowable F. 

H3. k^(T)>0, and k^CT) = 0 for i> N. No restriction is put on 

k^CT) for i = 2, ... ,N. 

THEOREM; C* furnishes the desired maximum of x-,(T). 

T T 



T 

Proof: We have shown that K^X 



= ^K^BF dt for every solution K to the 



0 0 



adjoint equation. Hence for the particular solution K and the matrix 
X* we have, on the path C*, 

(7) ki(T)x^(T) + k*(T)x2^ + ... + k*(T)Xj^^ = (K*^X*)^ 



, *T *s f ' 
= (K X )_ + \K 

° 0 






BF dt. 



Consider any other admissible path C' with F = F*. For this path and 
for K = K* we have 



(8) ki(T)x^(T) + ... + kjj(T)Xjj^ = (K ^X )q + ^BF dt. 

0 

Subtracting equation (8) from equation (7), we get 

T 

* * ' r '^’T t 

(9) kp(x^-Xj^) ^ UK BF -K BF ) dt ^ 0. 

0 

Hence, since k^(T)> 0, x^(T)^ x| at t = T. Hence the given hypotheses 
are indeed sufficient to give a maximim x^(T). Note that since kj(T) = 
0, for j=N+l, ... , n, the unspecified elements of X (T) play no part in 
the solution. 

Having developed a useful sufficiency condition, let us now consid- 
er the case where X is the matrix , where A and B are 2x2 matrices, 

\(diere F is the matrix and where we want to maximize x(T) subject to 

2 2 

the following constraints: First, (f]^) + (f2) = 1. Second, T>0 is 
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fixed. Third, solution curves must start at and end on the line y(T) 

y^, vThere y^ is given. For this problem the admissible curves are allow- 

able curves satisfying the third constraint. 

1 2 

Let us choose the solutions K and K to the adjoint equation having 
K^(T) =(^ and K^(T) = . Using these solutions in equation (6), we 



IT r IT 

(10) x(T) = (K X)q + \k BF dt 

0 

and 

T 

2T r 2T 

(11) y(T) = (K X)q + BF dt. 

0 

2 2 

Since the given constraint on F is that (f^) +(f 2 ) = 1, we may choose 

f^ = cos 0, f 2 = sin 0, where ^-^is a function of t. With this substitu- 
tion, a variation in 0 will cause a variation in F which will in turn 
cause a variation in x and in y. Let us call the variation in 9, 6 0 and 
the variations in x and in y at t=T, 6x and 6y, respectively. We then 
get, from equations (10) and (11), 

r IT 

(12) 6x = JK dt 

0 

and 

(13) 6y = ^k^^b 4|-60 dt, 

0 

for sufficiently small values of SQ* Let us now choose a special varia- 

IT B F 9T d F 

tion of the form 0 = e-^K B + ^2^ where e^ and e^ are con- 



stants yet to be chosen. Substituting this variation in equation (12), 
we get rjs 



(14) 6x = 



+ e2K^'^B-|-|-)dt 



12 



and a similar expression for 5 y. These can be simplified if we let 



i 



iT 

(K B- 



5 F 









-)dt. 



Using this substitution we get 

11 12 

(15) 6 X = e^I + e^I 

and 



(16) 6 y = e,I^^ + e„I 



r22 



12 21 22 
Note first that I = I . Note next the consequences if I =0. If 

92 2T ^ F 91 C 

I =0, then K B . ^ must be identically zero; hence I and by are 

both identically zero; hence y(T) is not affected by a change in the con- 
trol variable 6 . We describe such a situation by saying that the curve 
furnishes a stationary value for y(T). We will not want y to have a 

stationary value, and we shall henceforth assume for our curve that I 

9T ^ F 

0, or, equivalently, that B is not identically zero on the inter- 
val ( 0 ,t). 



Since T is fixed, equations (15) and (16) give the total differen- 
tials of X and of y. If we replace 6 y by y^ - y(T) , we can generally 
find values for e^ and which give a variation which will make the re- 
sulting curve admissible. 

Suppose now that we have an admissible curve and that we want to 
find conditions on the I^^ such that x(T) is a maximum. For an admissi- 
ble curve we have 

(17) 5 x = e^I + e 2 l 

and 

(18) 0 = + 62 !^^. 



Equations (17) and (18) can be solved for Sx>0 only if the rank of 
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6x 1 ^ 2 . I j 

is equal to the rank of 

0 l21 l22 ■ 



11 

21 



matrix is the same as the rank of 



'’T 

i22y 

/Sx 0 0 \ 

[ 0 l22j* 



But the rank of the first 



The rank of this last 



21 22 

matrix is one greater than the rank of the matrix (I I ); hence equa- 
tions (17) and (18) cannot be solved for 5x whenever the rank of 

^lll jl2\ 21 22 

I is equal to the rank of (I I ). But this will be true only 

l21 i22^ 

when the determinant j is zero. Hence if the admissible curve yields 

a maximum for x(T), then | | = 0* 

Let us consider some of the consequences of this condition. If the 
determinant = 0, then the matrix (I^*^*) has a zero eigenvalue. Let 



us choose constants c^ and C 2 such that 



(5) 



is an eigenvector of (I^*^) 



corresponding to this zero eigenvalue. Then let us consider 



(19) 



Clearly J^O. But J = + ( 02 )^!^^ 



= (c 



-P ■ 



■•■'P :»)(::)■ 






is an eigenvector corresponding to the zero eigenvalue of (I ), 

hence J = 0. Therefore c, ^ F + c^K^^B. ^ - — = 0. This equation 

1^0 2^0 

appears frequently in the literature on calculus of variations - it is 
known as the Euler equation . Curves whereon the Euler equation is sat- 
isfied are frequently called extremals . Hence K = c^K + C 2 K is a 



■^This is the definition given by Bolza 
given by Bliss [sj - 



Another definition is 
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solution to the adjoint equation which gives an extremal. 

•A* 

We saw that if c^> 0 and if 0 = 0 maximizes K BF then xCT) is a 

^ F 

maximum. The condition K B ^ ^ = 0 is a necessary condition for x(T) 

to be a maximum. We can and will choose c^> 0. We will also see that 

it is necessary that 0 maximize K BF; this is the Weierstrass condition. 
In this problem the Weierstrass condition is necessary and sufficient 
for an admissible curve to furnish the desired maximum. 

Let us assume that the rank of (I^^) is at least one for all curves 
which arise: this is called normal ity for the problem. Suppose now that 
we have an admissible, normal curve which furnishes the desired maximum. 
On it the Euler equation must be satisfied, and the matrix (I^^) has a 
zero eigenvalue with a corresponding eigenvector ; we will choose 

c^> 0. Then the F determined by the Euler equation maximizes the pro- 

duct K BF, as a function of 0, over the entire interval from t = 0 to 

t=T. The proof is as follows: 

Suppose that 0^(t) is the argument of F for some admissible curve 
satisfying the Euler equation but that 9^(t) does not maximize the pro- 
duct K ^BF(9) on some subinterval (t^,t 2 ) of (0,T). Then there is some 
other 9, say 92(t), such that K*^BF(92)> K*^BF(9j^) on the interval (t^, 
t 2 ). Let 9 = 92 on the interval (t^jt^ + dt^) , for t]^ + dt^<t 2 , and 
let 6 9 = e|K^^B-ll— + e 2 K^^B ^ ^ elsewhere on the interval (0,T) . 

Then, since both paths are admissible, 



(20) dx = K^^b(f(92) - F(9i)J dt^^ + + I^^e 2 

and 

(21) 0 = K^'^bCf( 92) - F(9;^)J dt^ + I^^e^ + 

If we now multiply equation (20) by Cj and (21) by C 2 and add, remember- 
ing that K*^ * ^1^^^ ^ 
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( 22 ) 



Cj^dx = K*^B £f(02^ “ f(0i^ dt^ 

since|^^jis an eigenvector of (I^^) corresponding to its zero eigen- 
value. Since the right-hand side of equation (22) and the constant c^ 
are both positive, the dx of equation (22) is positive. Furthermore it 
is possible to satisfy equations (20) and (21) by a set of e^^ which will 
give a positive dx and at the same time keep dy = 0. To see this, mul- 
tiply equation (20) by c^ and call the new equation (20 ); multiply 
equation (21) by C 2 and call the new equation (21). Add equation (20') 

to equation (21*) and call the resulting equation (23). Looking now at 

*T IT 2T 

equations (23) and (20), remembering that K = c^^K +< 02 ^ ’ have 

(23) Cj^dx = K*^B [^F(02) - F(0i)J dtj 

and 

(20) 0 = K^^B [f(02) - F(9^)J dt^ + I^^e^ + I^^e2. 

But, by assumption, K ^B [f( 02) - F(0i^>O, and / 0, hence we can 

solve these two equations for e^ and Hence it is possible to satis- 

fy equations (20) and (21) by a set of e^^ which will give a variation 
which will in turn yield a larger value of x(T). Hence if 9(t) is such 
that the product K*"^BF(0) is not a maximum on every subinterval of 
then x(T) will not be a maximum. It is wrthy of note that the form of 
the conditions for an extremum is independent of the particular problem. 

The Weierstrass-Erdmann corner condition is an immediate consequence 
of the above condition. For suppose that we have a curve whereon the 
hypotheses of the above theorem are satisfied. Suppose further that ti 

is a point at which the control variable §(t) is discontinuous. Then 

1 ^1 + 

^ =0. Proof: This must be so; otherwise the 

condition just found wuld not be satisfied in some neighborhood of t^, 
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K*’^BF 



*T I 

for either t> t, or t<t • The condition that K BF(0) = 

i i L 

(0) |x. is commonly known as the We iers trass -Erdmann Corner Condition. 
^1 + 

In this section we have introduced several important concepts as 
they pertain to a linear problem discussed in the following Sections. 

In the next section we consider a combination of linear systems, 
the combination being non-linear. 



V 
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III. The Weierstrass-Erdmann corner condition for a more general problem. 

In this section we consider a problem having a corner at a variable 
time. We develop the Weierstrass-Erdmann corner condition for this par- 
ticular problem. 

Let us consider the dif f erenliial equation 

(24) X = a'x + b’f 

where A’ is A, for 0< t< t^, and C, for t^< t<T, where B* is B, for 

0<t<ti9 and D, for t^< t<T X is the matrix^); F is the matrix 

A,B,C, and D are 2x2 matrices of functions of t which are piece- 

wise continuous and bounded on their respective intervals; B and D are 

non-singular. T is fixed; t^ is a variable to be determined. The con- 

9 2 

straint on F is as before, namely (f^) + (f 2 ) =1. In addition, x 

and y must be continuous at t^* For this problem the admissible curves 
are the allowable curves starting at X^ and ending on the line y(T) = y^. 

Note that this problem is non-linear on the interval (0,T), since 
A* and B* are functions of t^ as well as of t. The problem is, however, 
linear on each of the two sub-intervals (0,t^) and (t^,T). Hence we call 
this a semi-linear problem. 

We can rewrite equation (24) as 

AX + BF, 0 < t< ti 

(25) ^ 

_CX + DF, ti< t< T. 

The adjoint equation is 



(26) 



K = 



0< t< 

-c'^K, tj<t<T. 



Ve choose, as particular solutions to K = -C^K, tfe« solutions K^(T) = (g) 

9 /o\ 

and K'^(T) “Kj* For t^< t<T, then, each K is a function of T and of t. 
Since T is to remain fixed, however, we shall suppress it wherever it 
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occurs in the and shall consider the as functions of t only, on the 
interval (t^,T). Let us define a set of particular solutions to the 



equation K = 


= A^K by = K^(t^^), for i = 1,2. For t on the inter- 


val 


then, the are functions of the variables t^ and t where. 



as above, we suppress 

Before continuing further let us make one substitution. When con- 
venient, we shall use E to mean B, in the interval (0,t^) and D, in the 
interval (t^,T). (Note that E and B* are the same matrix). Using this 



convention. 


we get 


(27) 


T 

x(T) = K^^(tj^,0)X(0) + dt 

0 


and 


T 


(28) 


y(T) = K^^(ti,0)X(0) + ^K^^EF dt. 

0 



2 2 

Since the constraint on F is that (f^) + (f 2 ) = 1, we can replace f^ 

by cos 0 and f 2 by sin 0, where 0 is a function of t. Since we are con- 
sidering both 0 and t^ as variable, we get, from equations (27) and (28), 



(29) 


dx(T) = dt, +6x(T) 

d tj^ ^ 



and 



(30) 


dy(T) - dt, + 6y(T) 

o t^ 



where 8x(T) and Sy(T) are the variations in x(T) and in y(T) due to 
variations in the control variable 0. 

&x(T) and 6y(T) are given by the following equations for sufficient- 
ly small values of 0: 



(31) 


T 

6x(T) =^K^'^E-||- 60 dt 


and 


0 
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(32) 



66 dt 



6 y(T) = \k2Te-||- 






IT ^ F 2T d F 

Let us choose the special variation 59 = e^K ^ where 

e^^ and e,2 constants which remain to be specified. Then equation (31) 
becomes 

(33) 6x(T) = ^(K^'’^E-|-|-)(eiK^’^E-|-|- + e 2 K^^B|-|-) dt 



with a similar expression for 6y(T) . Adopting the substitution 



(34) I J ( K^"^E.|^) dt , 

0 

we get 

(35) 6x(T) 2 

and 

21 22 

(36) 6y(T) = e^I + e 2 l . 

We will use equations (35) and (36) later. 



The variations in x(T) and in y(T) due to a variation in t^ are 
more complicated. Rewriting equations (27) and (28) to show the way in 
which t^ enters, we get 



ti T 

(37) x(T) = K^^(ti ,0)X(0) + (K^^(t^,t)BFdt + ^K^^(t)DF dt 

6 ti 



and 



ti T 

(38) y(T) = K^^(tj,0)X(0) + ^K^^(t^ ,t) BFdt 



K^^(t>9P dt.. 



Hence 



(39) 



d t ^ 




d t^ 



(t^, 0 )x( 0 )dt^ 



+ 




(t^,t)BFdt^dt 



+ [(K^'^BF). -(K^^DF). J dt 

^ ^ 1 - ^ 1 +*^ 1 



and a similar expression for , the only difference being that 

d ti 
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IT 2T 

each K is replaced by K . To avoid problems in notation let us denote 

IT IT 

by 6K the variation in K ^ caused by varying t^ by a small amount dt^, 

It ^ 

i.e. to a first-order approximation OK = dt^^. With this nota- 

o t ^ 

tion, equation (39) becomes 

(40) dt^ = CK^^(tj,0)X(0) + y dt 



[(k1''bf)j^_-(k"df)^J dtj 



Since is not a function of t^ on the interval (t^,T), 6 = 0 on that 

interval. For t on the interval (0,t^), 6K^ satisfies the adjoint equa- 
tion, and hence 6K^ • = -A% on (0,tp. Keeping this in mind, let us 
consider the integral in equation (40). From equation (25), BF = X - AX. 
With this substitution, the integral in equation (40) becomes 



(41) 






6K^^(t^,t)X dt 






)AX dt. 



0 0 
Integrating the first integral by parts, we get 



IT 

(42) 6K (t^,t)X(t) 



- (6k . 



6K^ii>X dt. 



But 6K satisfies the adjoint equation, hence the integral in (42) is 

IT IT 

zero. Hence (42) reduces to 6K (t^,t^)X(t^) - (t^,0)X(0). But the 

IT 

variation in K due to a variation in t^, evaluated at t^, is equal to 
Qc^^(tj_)A - K^^(t^^)cJ dt^, by equation (25). Furthermore, X(t) is con- 
tinuous, so that X(tj ) = X(t]^+). Hence (42) is equal to 

[K^^(tj_)A - K^^(t^^)c3 X(t^) dt^ - 6K^'^(t^,0)X(0) , and equation (40) 
reduces to 



(43) 



Sx(T) 

d tn 



dtj = K^^(t^_) [AX(t^) + BF(t^_)J dti 
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- [cX(t^) + DF(t^^)J dt^. 



But, by equation (25), this sSys that 



(44) 



dt^ = [K^^(ti_)X(t^_) - K^^(t;^^)X(t^^)] dt^. 

d ti 

Similarly, 

(45) = [K^^(ti_)X(t^_) - K^^(t^^)X(t^^)J dt^. 

Hence equations (29) and (30) give, as the total variations in x(T) and 
y(T), ^ 

(46) dx(T) = ^ 60 dt - (k^'^x) | dt^ 

0 ^1- 



T ^ 

(47) dy(T) = { 60 dt - (K^'^X) I dt^ 

0 ‘^ 1 - 

IT 2T ^ F 

Choose the special variation 60 = (e-iK + e^ ) E . Then, using 

J- ^ 0 

equations (35) and (36), we get 



(48) dx(T) = I^^e^ + ~ 






dt. 



■ 1 - 



(49) dy(T) = - (K^'^X) 1 dt 



Observe that if t^ is fixed we have the problem discussed earlier. 

Let us assume, then, that admissible curves exist and that our prob- 
lem has a solution. If admissible curves exist, it is possible to find 
e^, and dt^ from equations (48) and (49) and hence to get a variation 
£0 such that the curve obtained by replacing 0 by 0 + 50 is admissible. 
Suppose that we have performed these calculations and have obtained an 
admissible curve. For this curve, dy(T) in equation (49) is zero. Be- 
fore continuing further, let us assume that I ^ 0. This assumption 
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is enough to insure normality in this problem. Having, then, an admis- 
sible normal curve, we can use the results obtained after equation (19). 
We showed that on an admissible normal curve for which x (T) is a 



maximum, the matrix (I^^) has one zero eigenvalue and that the components 
C]^ and C 2 of the eigenvector can be chosen so that C]^ is positive. We 
also showed that one solution to the adjoint equation which gives an ad- 
missible extremal is K = c^K + C 2 K"^. Hence let us choose c^ and C 2 as 
components of an eigenvector of the matrix (I^J) corresponding to the 
zero eigenvalue, with c^^^ 0. Multiplying equations (48) and (49) by 
these C£ and adding, remembering that the are continuous at t^, gives 



(50) Cjdx = K*'^(t 



j^X(t) dt^ . 



This leads to the following theorem: 

THEOREM : For x(T) to be a maximum, the quantity K^*^(t^) j^X(t) 

must equal zero. 

Proof: If the above quantity is positive, any dt^> 0 will give a posi- 
tive dx(T) and hence a larger x(T). If the above quantity is negative, 
any dt^<0 will give a positive dx(T). This is the Weierstrass-Erdmann 
corner condition for this problem. 

The next section uses the problem we have been studying for back- 
ground and considers the essentially different problem obtained by intro- 
ducing the constraint that x(t^) has some fixed value. We develop a 
numerical routine for determining the solution by a method of successive 
approximations . 
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IV. A numerical routine for determining the maximum x(T) for a fixed 
value of x(t^) 

In this section we consider the problem states as before but with 
the different constraint that x(t^) has some fixed value. We make up a 
curve consisting of two arcs, each of which is an extremal. We use the 
method of variation of extremals on these arcs to drive the resulting 
curve to admissibility and in a gradient technique to determine the curve 
on which x(T) is a maximum. 

Let us again consider equation (24), namely 
(24) X = a’x + B’F 

with the same conditions as in the beginning of Section III but with a 
different constraint on t^, namely that x(t^) = x^, where x^ is given. 

This problem, like the one in the preceding section, is semilinear. 
It is not linear on the entire interval (0,T), since a' and B* are 
functions of t^ as well as of t; it is, however, linear on each of the 
two sub-intervals (0,t]^) and (t^,T). 

The differential equation we have been working with is 



(25) 



X 



AX + BF, 0< t < t, 



\j:x + DF, tj< t< T' 

For this problem the admissible curves are allowable curves satisfying 
the above differential equation and such that X (|) = X^, x(tp = 
and y(T) = y^, where x^ and y^ are given constants. The adjoint equation 
for equation (25) is 



(26) 



K = 



j-A^K, 0 < t < tj^ 
tj< t< T. 



. 19 1 

We choose solutions K and K to the adjoint equation such that K (T) = 



0 



and K^'CT) 



0 



and such that ate continuous at t = t^. We further 



define the 2x2 matrix as before; its first column is K^(t) and its 
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second column is K^(t). Note that yc is a function of the variable t for 
t^<t< T and of both the variables t^ and t for 0<t<t^. 

Multiplying equation (25) on the left by , integrating from t = 0 
to t = t^, and using the fact that Tf is a solution to the adjoint equa- 
tion leads to 



(51) 



7^(t^,t)X(t) 




S ]f(t^,t)BF dt. 
0 



We have also, from the last section, 



fl T 

(37) x(T) = K^^(t^,0)X(0) + ) K^^(t^,t)BFdt + ^K^^(t)DFdt 



0 

^1 



^1 

T 



(38) y(T) = K^(t^,0)X(0) + ^ K^^(t^,t)BFdt + ^K^^(t)DFdt 



The differential equation is linear on each of the intervals. 

We choose F by using the method of variation of extremals on (0,t^) 
and each of them (t^,T). On each of the two arcs, we have the following 
by the theorem of section II. If C is an admissible arc, if K is the 
solution to the adjoint equation defined by K = c^K + C 2 K > where 
and are the solutions to the adjoint equation defined above, if c^> 0, 
and if, on C , F maximizes the scalar K EF, then C furnishes a max- 
imum X at the end of the arc, relative to the point at which the arc 
began. If we, then, choose F to maximize the product K ^BF, for 0<t 
<ti, and K DF, for t^ <t <T, the resulting curve will be made up of an 
extremal from t = 0 to t = t^ and an extremal from t = t^ to t = T. 

Once we have an admissible curve made up of two extremal arcs, we 
want some way to vary these curves so as to get a maximum x(T). We use 
a gradient technique, explained later, to choose a set of variations in 
the curve parameters to get a larger x(T) . The new curve may not be ad- 
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missible because of second-order effects. Hence, we again drive the 
curve to admissibility, say, by the method of variation of extremals 
mentioned above. ^ We continue this process until the maximum x(T) is 
obtained. 

The calculations proceed as follows: Since c^ is positive on each 

”^12 

arc, we can choose c^ = 1. Hence K = K + cK , where c is a constant. 
Unfortunately when a condition on X is given at time t^, may not be 
continuous at t^. Hence let us suppose that c is not the same in both 
arcs, and let us adopt the following notation: K = K + aK , for t on 
(0,t^), K* = + bK^, for t on (t^,T), where a, b, and t^ are unknowns 

which must be determined so that the curve made up of the arcs is admis- 
sible and yields a maximum x(T). Now consider F*. F* maximized K^BF 
on the interval (0,t^); on that interval K is a function of t^,t, and 
a. On the interval (t]^,T), K is a function of b and t only, hence F* 
is a function of b and t on the second interval. 

•k 

With F replaced by F , equations (51), (37), and (38) become 



(52) 

(53) 

(54) 



1^1 ^1 

■;t'^(t^,t)X(T) = 5 (ti,t)BF*dt 

'o 0 

ti T 

x(T) = K^^(t]^,0)X(0) + ^K^’’^(t^,t)BF*dt + 



)DF*dt 



2T 



y(T) = K (t^,0)X(0) 



ti T 

+ ^K^^(t^,t)BF*dt + ^K^’’‘^(t)DF*dt 



From these equations we want to devise a routine that will, first, give 



^It is possible that we might obtain an x (T) from a set of curves 
that are admissible under one admissibility criterion but not under 
another. After trying various admissibility criteria, the author de- 
cided upon that of calling a curve admissible if 
$ 10 "^. 



lx(tj)-x^| + |y(T)-yf\ 
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US a set of admissible curves and, second, find an admissible curve where 
on x(T) has its maximum value. We want also to find a condition which 
will indicate that no admissible curves exist, if such is indeed the case 
In this problem, a curve can be characterized by a set of values for t^? 
a, and b* Hence we want a routine to find values of a, b, and t^ which 
will determine an admissible curve whereon x(T) is a maximum. 

Let us first get expressions for the total differentials of x(t^), 
of x(T), and of y(T). First, of x(t^) . Equation (52) can be written 

ti 

(55) = 7C^(t^,0)X(0) + J^(t^,t)BF*dt. 

0 

Taking differentials of both sides gives 

(56) d;^'^(t^,t^)X(t^) + ;^'^(t^,t^)dX(t^) = d-^'^(t^,0)X(0) + 

tl 

:^^(ti,t^)BF*(t^_)dt^ + J d ^'^(t^,t)BF*dt + 

0 

rl T* * 

) (t^,t)B6F dt. 

To get d (t^,t^), note that the first t^ denotes the end of the inter- 
val (0,t^) and that the second is the value that the running variable t 
assumed at the end of the interval. The differential due to the change 
in the first t^ is, since satisfies the adjoint equation, ^^(t^,t^) 
(A-C)dt^; the differential due to the change in the second t^ is -l^^(t^, 
t^)A dt^. Hence d^^(t^,t^) = - )^(t^ ,t^)Cdt and d ^^^(t^ ,t ^)X(t^) = 

— ^'^(t 2 , t^)CX(t^)dt^ . 

T * 

Consider next V d 3^ (t^,t)BF dt. By the argument after equation 
0 

(42), this integral is equal to A^(t^ ,t^) (A-C) X(t^)dt^ - d^^(t^ ,0) X(0) . 
Substituting the above results into equation (56) yields 
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(57) '^'^(t^,t^)dX(t^) = ■^'’^(t^,t^)AX(t^)dt^ + 

;^:'^(tpt^)BF(t^_)dt^ + J /[:'^(t^,t)B6F*dt 

0 

= ^^(t^,t^)X(t^_)dt^ ^^(t^,t)B6F*dt. 

0 

Calculating the integral in the above equation is somewhat more diffi- 

■k IT 2T 

cult* F is chosen to maximize the produce (K + aK )BF. Expanding 

IT 2T 

the product (K + aK )B yields the matrix 

[kll»^12Jbj2-t''*ak22jb2) . 



Calling the first element h^ and the second h^ and remembering both that 
F must maximize the product (h^ h 2 >F and that F can be written as 

/cos 9\ that cos 0 = h-^/ + (h 2 )^J ^ and that sin 6 = h 2 / 

\sin 0/ ^ . 

r 2 01 i< * /-sin 0 \ -1 

[jCh]^) + C^ 2 ^J Hence 6 F = f 0/^^’ ® (.h^/h^) . 

Substituting in the expressions for sin 0 and for cos 0 and finding 60 
in terms of 6 h]^ and 6 h 2 yields 



* 



6 F 



f-h2\ h;^8h2 - h2&h^ 

Ul/ '[(hp^ .(h2)^J^^^ ’ 



6 h^ and 6 h 2 are each the result of variations both in t^ and in a. Con- 
sider first the part of the variation (h^ bh^-h^^h^) due to a variation 
in a. 

(58) = [(k“*.k^2)b„ * <k21*ak22)b^;[fcl2bj2 . k^%3 

- [(k^^+ak^^)b^2 * 

+k^2b2i)3 + j(k^H^^+k21b2i)(k^^b^2-"'^^^^22^ " ^*^^^^12 
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I 



1 21, . ,, 12, , 22, n 

+k b22)(k + k b^^)! 

. 11 , 12 ,^ , u V N 1 21 , 22 K V S 

— k k X 1 ^12~^12^1 k k ^^21^22~^22^21^ 

+k^V^(b^^b22-b^2^21^ ^ k^^k^^(b2ib^2-‘’22^lP 

= (k k -k k ) (b^ 3 ^b 22 -b^ 2 ^ 21 ^ 

= |/^| • | b | , 

where 1 7^1 denotes the determinant of ^ and similarly for jsl . Hence 
the variation due to a variation in a isj/^| |b| da. 

The part of the variation (h]^6h2-h2&hj^) due to a variation in t^ is 

3 h 2 3 h^ p ^2 21 22 “i 

( 59 ) - ^ 2 ^^ *^^1 = +ak^^)b^^+(k^ +ak^^)b 2 J 

*[( 6 k^^ + a 6 k^^)b^ 2 + (&k^^+a 6 k^^)b 22 ]- 
fCk^^ + ak^^)b^2 + (k^^+ak^^)b22 j {(6k^^ 
+a£fc^^)b,, + ( 6 k^^ + a 6 k^^)b 3 



= B 



'11 

.11 .12 :;, 11 „ 12 

k + ak : rfk + a6k 

21 22 / 21 22 

k^”^+ ak ^k + a6k 



For convenience, we shall refer to this variation as VAR(l) henceforth. 

Note that each of the bk^*^ can be calculated, since 

# 

7^(ti,t^)dt^ and since 6 J^= -A^6/S^on (0,t^). 

T 

It is possible to simplify, within the integral, the terms Al^(tj,t)B 
^”1^2 j, as follows: 



7 (f(ti,t)B ' 




,11 21 \ / , \ / ,, 11 , 12 . , ,, 21 , 22 _ 
k k \ /bj^Q^ *^12 \ l~ )b^2'*'(^ )b 



k ^2 U21 b 22 



,^ 2 -v^ ^..22 

(k^^+ak^^)b^^+(k^^+ak^^)b 2 i 



= (■;) H 14 
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Using the above substitutions, equation (57) simplifies to 



(60) -^T(t^,ti)dX(t^) = -^’^(t^,t^)X(t^_) + 



r 

\ HiBiiKi H 

\ 1/ r<h^)2 + (h2)2j 3/2 



dt. 



is a non-singular square matrix, hence ^ exists, and 

(61) dX(t^) = X(t^_) dt^ + 



-1 



ailelM 



|b| da + VAR(l) dt 



1- dt, 



where dX(tj^) 



as 



K )2 » 

/dx(t^)\ ^ future reference, let us rewrite dx(t,) 

Vdy(ti)y 1 



(62) dx(t^) = + Oi^ 2 ^ay 

^ x(t]^) d X(t]^) 

where , = — and a, « = — r . Thus we have derived an ex- 

11 <5 t^ 12 a 

pression for dx(t^) in terms of variables which we can calculate. 

We want next to derive an explicit expression for dx(T) . We can 
rewrite equation (53) as 



T 

(63) x(T) = K^'^(t^,tj)X(t^) + 5 K^’’^(t)DF*dt, 



whence 

IT IT 

(64) dx(T) = dK^\t^,t^)X(t^) + K (t^,t^)dX(t^) 

T 

IT * /* IT * 

- K (t)DF (t^^)dt^ J ^ (t)D6F dt. 

IT IT 

By the discussion preceding equation (57), dK (tj^,t^) = -K (t^,t^)Cdt^, 
and dK^^(tj^ ,t^)X(t^) = -K^'^(t , t^)CX(t ^^)dtj^ . If we write the differen- 
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tial for dy(t^) corresponding to equation (62) as dy(tj^) = a^^dt^^ + 
the product K^^(t^,t j^)dX(t^) becomes 

By an argument similar to that preceding equation (60), noting that F 

is dependent only on b in the inteirval (t]^,T), we get 
T T 

(65) \ K^'^(t)D6F*dt = \ (-b) IbI ^ jtCpd tidb 

£hi)2+(h2)2]^/2 

Combining the above results yields 



(66) dx(T) = -K^^(t^jt^) ^X(t^)+DF (t^^)Jdt^ +(k^ V^^+k^^a^^)dt^ 

T 

+ (k^^ 0 L„+k^^a )da+ \(-b) 1 b| ^ (fCpdt db- 

Combining terms gives us an equation corresponding to equation (62), 
namely 

(67) dx(T) = 0^21^^! where^^i ^ , 

"22 = - 1 ^ ■ "32 = . 

dy(T) can be derived much as we derived dx(T); the equation corresponding 
to equation (66) is 

2T, . r- . * . 



11^ , 12. 



( 68 ) 



dy(T) = -K^^(t^,t^) [cX(tj) + DF (t^^)^dt^ 



+ (k^\i^^+k^^agj^)dt + (k^ 'a , ^+k^^a^^)da 



+k22a 
12 & 2 ' 



+ f (d(^ j/CP d t db, 

J[(hp2,(hp2j3/2 



which we abbreviate as 



(69) dy(T) =» 3 idt^ + a^^'^a +® 33 db. 
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We thus have three equations to use to determine corrections in t^ 
and b which will lead to a maximum x(T), namely 



(62) 


dx(t j) 






(67) 


dx(T) 




+ a da + a 

22 


(69) 


dy(T) 

J 


= “3l‘^S 
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The numerical procedure used for solving for a maximum x(T) is the 
following: Choose values for t^, a, and b. With these values, calcu- 

late xCtj), x(T), and y(T) and theol^^, starting from X(0) = X^. Prob- 
ably the curve so determined will not be admissible, i.e. x(t^) will not 
be x^ and y(T) will not be y^. To obtain a curve which will be admissi- 
ble, solve equations (62) and (69) for two of the variables dt^, da, or 
db, setting the third equal to zero. It is possible to solve equations 
(62) and (69) for dt^, da, and db only if the rank of the matrix 



(70) 



'a 



11 “12 



“31 ®32 “ 33 ' 



is two, hence this is a criterion for deciding whether any admissible 
curves exist. If the rank of (70) is two, replace a, b, and t^ by a + da, 
b + db, and t^ + dt^ and repeat the calculations. Continue this process 
until either the rank of (70) becomes one or until some admissibility 
criterion is met. Note: In correcting two variables, it is essential 

that the value of third variable be chosen close to the value it will 
have on an extremal; otherwise the correction routine may not converge 
to an admissible path, even when such a path exists. The admissibility 
criterion used by the author was |x^ - x(t 2 ^)j + jy^ - y(T)| f 10"^; for 

the matrices used this criterion was normally met in less than ten iter- 
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at ions when the determinants 



‘^ll ^^12 


and 


‘^ii 0 


°-31 “32 




°^31 ®33 



went to zero, this 



was also obvious by the tenth iteration. 

When admissibility had been achieved, equations (62), (67), and 
(69) took on the form 

(71) 0 = a^^dt^ + 



0 = + 01^2^^ + a^^db. 



(72) dx(T) = a 2 i^t^ + a 22 ^^ ^23^^ 

(7 3) 

These equations were solved for dx(T) unless the rank of 

ail «i2 ® \ 



tt2i 022 “23 



'^1 “12 ° \ 
^®31 “32 “33) 



was equal to the rank of [ - - \ . If the rank of 

S^^31 ®32 ^33/ 

the latter matrix was two, this condition was that the determinant of the 
first matrix must be zero. We therefore want some way to determine a new 
set of corrections to drive the determinant toward zero. The method 
adopted was to choose dt^, da, and db proportional to the components of 
grad x(t^) X grad y(T) , where indicates the vector cross-product and 
where the proportionality constant is some fixed number e multiplied by 
the determinant j | • The old values for and b were replaced by 

the corrected values and a new admissible curve was found using the meth- 
od outlined above. For this curve the determinant correc- 

tions da, db, and dt^ were again calculated together with x (T) . This 
process of finding an admissible curve, then calculating a set of cor- 
rections to increase x(T), then finding a new admissible curve, then 
finding a new set of corrections to increase x(T) was continued until 
the determinant changed sign. When the determinant changed sign, 

however, the above procedure no longer worked - the corrections became 
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larger rather than smaller, and the values of x(T) got smaller instead of 
larger. A procedure that worked when the determinant changed sign was 
the following: A new proportionality constant was chosen equal to one 

fifth of the sum of the positive and negative determinants; a new set of 
corrections da, db, and dt^ was then calculated. The admissible curve 
calculated after finding this new set of corrections invariably yielded 
a larger x(T) than had been obtained before. This method of successive 
approximations seems to be a nearly foolproof method for determining the 
maximum x(T) . The procedure was stopped when either the absolute value 
of the determinant was less than 10"^ or it became obvious that the con- 
vergence was too slow for the amount of time available on the computer. 

If convergence was too slow, it could usually be improved by increasing 
the number £ referred to above. 

Another successive approximation procedure used was to compare the 
x(T) resulting from replacing t^,a, and b by their corrected values with 
the previous x(T). If the new x(T) were smaller, £ was reduced by a fac- 
tor of ten and a new set of corrections and hence a new x(T) computed. 
This process was continued until either the sum of the absolute values 
of the corrections was less than 10“^ or the new x(T) was larger than 
the old x(T), in which case the new x(T) became the value compared with. 
This process gives a monotonically increasing sequence of values of x(T); 
the speed of computation was comparable to that of the other method. 

It was necessary to insure that the corrections from the routine 
were not so large as to make the linearity approximations in the correc- 
tion integrals invalid. The author used the following procedure: In 

the correction to admissibility, if either da or db exceeded .3 in abso- 
lute value, or if dt^ exceeded .5 in absolute value, he divided all the 
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corrections in half. For the problems worked, this criterion was ade- 
quate although in several cases it slowed convergence quite a bit. (For 
example, where the initial value given for a was 1. and where the final 
value on the admissible curve was 12.8.) Once admissibility was achieved, 
it was necessary to use some linearity criterion on the second set of 
corrections. The procedure was to divide the corrections in half if the 
sum of the absolute values of the corrections exceeded .5. 

With the above correction procedures, computations were made on 
several sets of matrices with various boundary conditions. The only case 
where it was possible to verify the results by comparison with a physical 
situation where A and C were set equal to zero and where B and D were 

2x2 scalar matrices. This set of matrices gives the differential equa- 

faF, 0 <t <ti 9 0 

tion X = 7 „ ^ m • constraint on F was that (fn) + (f^)^ = 1; 

(jc. F,t^<t<T ^ 

f^ and f 2 were consequently chosen as cos a and sin a, for 0 < t< t^, and 

as cos p and sin p, for t^<t<T. Solutions to the differential equation 

as set up describe the path taken by a light ray going from one isotropic 

medium into another, with a and P being the angles between the light ray 

and the normals to the plane of discontinuity. As such, these solutions 

f 

should obey Snell’s Law, namely that the ratio of the velocities on oppo- 
site sides of the discontinuity plae should be the same as the ratio of 
the sines of the angles the ray makes with the normal to that plane. The 



cases tried were for = l.Oy y^ = 5.0, a = 1.0 and b = 2.0 

and, in the second case, a = 2.0 and b = 1.0. For the first case, the 
ratio between the sines should be .5. a was found to be approximately 
.403378 radians and ^ approximately .902758 radians. The<*'ratio of the 
sines of these angles was .500, which verified the accuracy of the rou- 
tine. For a = 2.0 and b = 1.0, the angles were reversed, and again the 
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accuracy of the routine was verified. 
Another case tried was where A = 

/q o\ /2.0 0. \ 

= (0 0)’ ° = (0. 2.oj’ 



' 1.0 ,l\ _/l.O .2\ 

^1.2 l.oj’ ® .0 1.1/ ’ 

= 1.0, T = 4.0, and where initial 



values for t^, tana? and tan p were .73, 4.0, and .7, respectively. 

The routine converged to t^ = .8935, tan a = 1.4925, and tan P = .4967, 
with the maximum x(T) being 6.4622. This convergence took fifteen min- 
utes and fifty-one seconds on a GDC 1604 computer. The author felt that 
convergence could be improved if better initial estimates were made of 
tan a, and tan p. 
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appendix I 




These constants are given: 

x(0), y(0), XMID. YHNAL. MATRICES A. 

B. C. and D. EPSFAC. XU. ATC12. ATDll 
TONE, COUESS. DGUESS. CAPT. ERROR. 
EPSl 




SET STEP SIZE 



TSTEP(l) * (CAPT - TONE)/50 
TSTEP(2) = TON E/50 



1 

INT^GRi^ION FROM^ - Q_APT to L*_TONE 
In: /r (CAPT) *(JJ) 

Integrate: K ■ -C^A' 

OuT:"^7t7) 



Calculate: 6 A' (tj) = (A*^ - c’^)A' (t^) 



^NT^GRATON L*-?_ 

In:/r(t^). 6/T(tp 

Integrate: ^ * -A^A', 6 A' = ~A^6K 



Out: K (0). 6A'(0) 



I 



INTEGRATION FROM t = 0 to t =* TONE 






-A^K, 6/r * -A^6K 



Integrate: A' 

Calculate: (aCON( 1) AC0N(2)) » + CGUESS x B 



SRT = ^(aCON(1))^ + ^ACON(2)^^^ 



1/2 



r - (Sg!) ^ 

Integrate: X = AX + BF 

d (cORR(l)^ = (-CGUESS x |bJ^ “ K| ^ ’ (SRT)®) dt 

d (;ORR(2 )) = (|b|® X |/r|® +(SRX)®)dt 



Calculate: OINTFAC - |B| 



NOTES ^ 

I.FLOVV CHART PAGE NUMBERS ARE 
DESIGNATED BY THE LETTER F, e.g., 

FG IS page S of the FLOW CHART. 

aCONNECTC.RS ARE REFERENCED BY 
PAGE NUMBERS, e.g..(T^ f2 DESIGNATES 
, CONNECTORS IS ON'^PAGE F2. j 



+ CGUESS X 



+ CGUESS X 6k“ 



+ CGUESS X K®® 6K®^ ♦ CGUESS x 6K®® 



Integrate: d (C0RR(3)) = -CGUESS x |b|® x |/r| x DINTFAC + (SRT)® dt 

d (cORR(4)) = )b|® X |A"| X DINTFAC + (SRT)® dt 



i 



Set XINMID* x(t^) 



Calculate CKI » 



CHANGES IN X(t^) AND Y DUE 
TO CHANGES IN CGUESS 
XONE(2) - (10)CKl(Cg“U)) 
YONE(2) - CORR(2) 



T 



VARIATION IN X(TONE) DUE TO A VARIATION > 
IN TONE 

XONE(l) - (1 0) CKI (cOTR(4)) * XDOTI 



T 



/PARTIAL VARIATION IN Y DUE\ 
TO A VARIATION IN TONE 



YONE(3) 



CORR(4) 



R: X(t^J - (ytoTi)* CORR(l). CORR(2) . CORR(3) . CORR(4) 



O" 



38(F)) 





39 (F2) 




40 (F3) 






QdeterI 



1 



s 10'® ? 



3 



No 






Solve for DELTAD and DELTAT 

/xdiffN /xONE(1) 0 \ 

^YDIFF^ “ \V0NE(1) CORA ^ 


/deltat\ 

\DELTAD/ 


0- 


> 






(IdeltadI < 


.3 ? 








^No 





Solve for DELTAC and DELTAD 


fxDlFF] 

^ydiff/ 


. /X0NE(2) 
VC0RR(2) 


0 N /deltacX 
CORA ) \DELTADy 








XONE(2) X CORH(2)| s 10 ? 



5 « io 



Divide all corrections by 
two 



5 



F4 



Ves 



Print: NO ADMISSIBLE 
PATHS EXIST 




" > ^DELTACl <, 







Yes 



"“"V No 



Divide all correc* 
tions by two 



10 JF4 



± 



]deltat| < 5 ? 

I No 



F2 



Divide all corrections by 
two 






F4 
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APPENDIX II 



PROGRAM D I SCON 

This program computes max x(Capt) for a given y(capt> and a given 

X(T1), WHERE XDOT»AX+BF,FOR T BETWEEN 0 AND T1*AND CX+DF,FOR T 
BETWEEN T1 AND CAPT. CAPT IS GIVEN, BUT T1 IS NOT. 

DIMENSION A<2»2),B{2,2)»C(2»2),D(2,2)»BK(2»2),F(4),TSTEP(2),E<2,2, 

1 2)»ELV(14),DE(16), AAK ( 2 ♦ 2 ) ,DP ( 16 ) , ACON ( A ) . AKTB ( 2 , 2 ) ,CORR 

2 (4) ,CKI (2»2) »FLES< 2 ) ,XONE( 5) ,Y0NE<4) , B<A ( 2 , 2 ) , FMOR ( 2 ) ,DELK(2,2) 
3 »DACON( 4) ,DFE(4) ,AK( 5,16) 

READ 1503,EPSFAC 
1503 FORMATC 1F20.10) 

READ 501,((A(I,J),J»1,2),I=1,2)»((B(I,J),J*1,2),I=1,2) 

501 FORMAT (8F10.1) 

READ 501, ( (C( I , J) » J=l,2 ) » 1=1 »2) ♦ ( (D( I , J) ,J=1 ,2 ) , 1=1 »2 ) 

PRINT 551 

551 FORMAT ( 14X , IHA , 29X , IHB » 29X , iHC » 29X , IHD ) 

PRINT 552, ( (A( I ,J) ♦J = l »2 ) ♦ <B( I , J) ,J = 1,2 ) , <C( I , J) , J=l,2) ,( D( I ,J) ,J = 

• + 1 , 2 ), 1 = 1 , 2 ) 

55 2 FORMAT {4{8X,F4,1,6X,F4.1,8X) ) 

READ 502,TONE,XO,YO,YFINAL,XMID,CGUFSS,ERROR,CAPT 

502 FORMAT (8F10.8) 

READ 503,DGUESS 

503 FORMAT( F10.8) 

PRINT 553,TONE,XO,YO,YFINAL,XMID,CGUESS,ERROR,CAPT 
553 FORMAT! 4X,12HT0NE GUESS= ♦F5.2» 9H ( XO, YO ) = ( , F4 . 1 , 1H , , F4. 1 , 9H)YFIN 
+AL= ,F4.1»7H XMID= »F5.2,9H CGUESS* ,F4.1,8H ERROR* ,F7.6»6H CAPT= 
+,F5.2) 

C INTEGRATE ADJOINT BACKWARDS TO GET K(0) AND KINVERSE AT Tl, 



001 
OOt . 
00< " 
OOt 
OOt 
OOi 
not 
OOi 
00 < 
00 < 
on I 
00 ! 
00 
00 
00 
00 
00 
00 
00 I 
00 ! 
oo:, 
oo; 
00 ; 
oo; 
00 . 
oo;i 

oo: 



F(1)=0. 00. 

F(2)=0.5 00 
F(3)=0.5 oo; 
F(4)=l, 00' 



ATC11=5. oo: 

ATD11=5. oo: 
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I 

t 



Xll=l. 

EPS2=-.01 

1 TSTEP( 2 )=TONE/50. 

DO 1003 I=l»3 
YONE( I )=0. 

003 XONE(I)=0. 

TSTEP( 1)=(CAPT-TONE)/50. 

TSTEP( 1 )=-TSTEP( 1 ) 

TSTEP{2)=-TSTEP(2) 

BK(1,1)=1. 

BK(2,1 )=0. 

BK( 1.2)=0. 

BK(2»2)=1. 

T=CAPT 

DO 101 1=1,2 
DO 101 J=l,2 
E(2,I,J)=A( I ,J) 

t lOl E(1,I,J)=C{ I,J) 

ELV{1)=BK(1,1) 

ELV(2)=BK{ 1,2) 

ELV{3)=BK{2,1) 

ELV{A)=BK(2,2) 

DO lOOA IMA=5,8 
1004 ELV(IMA)=0. 

IB=1 . 

1001 DO 102 IE=1,50 
DO 103 IC=1,4 
DO 104 ID=1,4 

104 DE( ID)=ELV( ID)+F<IC)*AK(IC-1,TD) 

DP(1)=-E( IB,1,1)*DE< 1)-E( IB,2,1)*DE(3) 
DP(3)=-E( IB,1,2)*DE( 1 )-E( I B, 2 , 2 ) *DE ( 3 ) 
DP(2)=-E( I8,1,1)*DE( 2>“E( IB,2,1)*DE(4) 
DP(4)=-E( IB,1»2)*DE<2)-E( I B, 2 , 2 ) *DE ( 4 ) 
DO 103 ID=1,4 

103 AK( IC, ID)=TSTEP( IB)*DP( ID) 

^ DO 105 ID=1,4 



0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 

0053 

0054 

0055 

0056 

0057 

0058 

0059 

0060 
0061 
0062 

0063 

0064 

0065 

0066 

0067 

0068 
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10 5 ELV( I0)=ELV( I D ) + ( AK ( 1 . I D ) +2 . *AK ( 2 » I D ) +2 . *AK ( 3 . I D ) +AK ( 4 , I D ) ) /6. 0069 

102 T=T+TSTEP( IB) 0070 

BK(1»1)=ELV(1) 0071; 

BK( 1»2)=ELV(2) 0072 

BK(2»1)=ELV(3) 0073 

BK(2.2 )-ELV(4) 007« 

IB = 2 007? 

1002 ELV(5)=(A( 1»1)-C( 1.1) ) *BK ( 1 » 1 ) + ( A ( 2 » 1 ) -C ( 2 » 1 ) )*BK(2.1) 007f 

ELV(6)= (A( 1,1 )-C( 1,1) )#BK( 1 ,2)+( A(2,l )-C(2,l ) )*BK(2,2) 007- 

ELV(7)=(A(1,2)-C( 1,2) ) *BK ( 1 , 1 ) + ( A ( 2 , 2 ) -C ( 2 , 2 ) )*BK(2,1) 00 7f 

ELV(8)=(A(1,2)-C( 1,2) ) *BK ( 1 , 2 ) + ( A ( 2 , 2 ) -C ( 2 , 2 ) )*BK(2,2) 007? 

DO 112 IA=1,50 008f 

DO 903 IC = 1,4 0083i 

DO 904 ID = 1,8 008: 

904 DE( ID)=ELV( ID)+F( IC)»AK( IC-1,TD) ' 008: 

DP(1)=-E( IB,1,1)*DE(1)-E(IB,2,1)*DE(3) 008', 

DP(3)=-E( IB,1,2)*DE( 1)-E( I B , 2 , 2 ) *DE ( 3 ) 008' 

DP(2)=-E( IB,1 ,1 )*DE( 2)-E( I B , 2 , 1 ) *DE ( 4 ) 00 R< 

DP(4)=-E( IB,1,2)*DE(2)-E( I B , 2 , 2 ) *DE ( 4 ) 008' 

DP(5)=-E( IB,l,l)*DE(5)-E( I B , 2 , 1 ) *0E ( 7 ) 0081 

DP(7)=-E( IB,1,2)*DE(5)-E( I B , 2 , 2 ) *DE ( 7 ) 008‘ 

DP(6)=-E( IB,1,1 )*-DE(6)-E( I B , 2 , 1 ) *0E ( 8 ) 009(1 

DP(8)=-E( IB,1,2)*DE(6)-E( I B , 2 , 2 ) *DE ( 8 ) 009' 

DO 903 ID»1,8 009. 

903 AK( IC, ID)=TSTEP( IB)*DP( ID) 009 

DO 905 ID=1,8 009' 

90 5 ELV( ID)=ELV( ID) + (AK( 1 » I D) +2. *AK ( 2 , 1 0 ) +2 . *AK ( 3 , I D ) +AK ( 4 , I D ) ) /6. 00 9 

112 T=T+TSTEP( IB) 009 

BK( 1 ,1 ) =ELV( 1 ) 009 

BK( 1,2)=ELV(2) 009 

BK(2,1 )=ELV(3) 009 

BK(2,2)=ELV(4) 010 

DELK(1,1)=ELV(5) , 010 

0ELK(1,2)«ELV(6) ’ 010 

DELK(2,1)=ELV(7) ' 010 

DELK(2,2)=ELV(8) 010 
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DELK IS THE MATRIX OF DELK S AT T»0 OIOS 

TSTEP( 1)=-TSTEP(1) 0106 

TSTEP( 2)=-TSTEP(2) 0107 

THE K MATRIX, AS IS, IS AT T=0, 0108 

THE NEXT INTEGRATION IS FROM T»0 TO T=TONE 0109 

ELV(1)=BK(1,1) 0110 

ELV(2)=BK(1,2) 0111 

ELV(3)=BK(2,1) 0112 

ELV(4)=BK(2,2) 0113 

ELV(5)=0. 0114 

ELV(6)=0. 0115 

ELV(7)=X0 0116 

ELV(8)=Y0 0117 

ELV(9)=DELK(1,1) 0118 

ELV(10)=DELK(1,2) 0119 

ELVdl )=DELK(2,1) 0120 

ELV(12)=DELK(2,2) 0121 

ELV(13)=0. ^ 0122 

ELV(14)=0. ' 0123 

> DO 202 IA=1,50 0124 

DO 203 IC=1,4 ' 0125 

DO 204 ID=1,14 0126 

204 DE( ID)=ELV( ID)+F( IC)*AK(IC-1,ID) 0127 

DP(1)=-A(1,1)*DE(1)-A(2,1)*DE(3) 0128 

I DP(3)=-A(1,2)*DE(1)“A(2,2)*DE(3) 0129 

DP(2)=-A(1,1)*DE(2)-A(2,1)*DE(4) 0130 

0P(4)=-A(1,2)*DE(2)-A(2,2)*0E(4) 0131 

i AC0N(1)= DE(1)*B(1,1)+DE(3)*B(2,1)+C6UESS*(DE(2)*B( 1,1)+DE(4)*B(2, 0132 

' +1)) 0133 

AC0N(2 )=DE( 1 )*B( 1,2)+DE(3)*B(2,2)+CGUESS*(DE(2)*B(1,2)+DE(4)*B( 2, 0134 

I +2)) 0135 

' SRT=SQRTF(ACON( 1)*AC0N(1)+AC0N(2)*AC0N(2) ) 0136 

1 AKDET=DE(1)*DE(4)-DE(2)*DE(3) 0137 

BDET=B( 1,1)*B(2,2)-B(2,1)*B(1,2) 0138 

. DP(5)=-CGUESS*AKDET*AKDET*BDET*BDET/SRT*»3 0139 

' DP(6)* AKDET*AKDET*B0ET*BDET/SRT**3 0140 

I 
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DP(7)=A(1»1)*DE(7)+A( 1 , 2 ) *DE ( 8 ) + ( 6 ( 1 ♦ 1 ) *ACON ( 1 ) +B ( 1 ♦ 2 ) * ACON ( 2 ) ) /SR OlAl 
+T 0142 

DP(8)=A(2»1 )*DE(7)+A(2»2)*DE(8)+(B(2»1 )*AC0N( 1 ) +B { 2 ♦ 2 ) *AC0N ( 2 ) ) /SR 0143 
+T 0144 

DP(9)=-A( !♦ 1 )*DE(9)-A(2»1 )*DE( ll) 0145 

DP( 11 )=-A( 1 t2 )*DE( 9)-A (2.2 )*DE( 11 ) 0146 

DP(10)=-A( 1.1)*DE( 10)-A(2.1)*DE( 12) 0147 

DP( 12)=-A( 1 .2 )*DE( 10)-A(2»2)*DE( 12 ) 0148 

DINTFAC=CGUESS*(DE( 1)*DE(12)-DE(4)*DE(9)+DE( 2)*DE( 11)-DE(3)*DE( 10) 0149 

1 +CGUESS*(DE(2 )*0E(12)-DE(4)*DE( 10) ) ) -DE ( 1 ) *DE ( 1 1 ) +0E ( 3 ) *DE ( 9 ) 0150 

DP( 13)=-CGUESS*BDET*BDET*AKDET*OINTFAC/SRT**3 0151 

DP(14)= BDET*BDET*AKDET*DINTFAC/SRT**3 0152 

DO 203 ID=1 ,14 0153 

203 AK( IC, ID)=TSTEP(2)*DP( ID) 0154 

DO 206 10=1,14 0155 

206 ELV( ID) =ELV( I D ) + ( AK( 1 , I D ) +2 ,*AK ( 2 , I D ) +2 ,*AK ( 3 , I D ) +AK ( 4 , I D ) ) /6. 0156 

BK(1,1)=ELV(1) 0157 

BK( 1 ,2 )=ELV( 2 ) 0158 

BK(2,1 )=ELV( 3) 0159 

BK(2,2 )=ELV(4) 0160 

CORRd )=ELV(5) 0161 

C0RR(2 )=ELV(6) 0162 

XM0D=ELV(7) 0163 

YM0D=ELV(8) . 0164 

DELK(1,1)=ELV<9) 0165 

DELK(1,2)=ELV(10) 0166 

DELK(2,1 )=ELV( 11) ' 0167 

DELK(2,2)=ELV( 12) 0168 

CORR(3)=ELV(13) 0169 

C0RR(4) =ELV( 14) 0170 

202 T=T+TSTEP(2) 0171 

CORR 3 AND 4 GIVE THE VARIATIONS OF X AND Y AT T1 DUE TO DTI 0172 

GET K INVERSE AT TONE AND F MINUS 0173 

AKDET=BK(1,1)*BK(2,2)-BK( 1,2)*BK(2,1 ) 0174 

XD0T1=DP(7) 0175 

YD0T1=DP(8) 0176 
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I 



I 

L CKI (1.1)=BK(2»2)/AKDET 0177 

I CKI (1*2)=-BK(2,1)/AK0ET 0178 

CKI (2, 1 )=-BK( 1 »2 ) /AKDET 0179 

CKI <2*2)=BK( 1»1)/AKDET 0180 

i FLES(l)=ACON(l)/SRT 0181 

f FLES(2 )=AC0N(2)/SRT 0182 

Y0NE(3)=C0RR(4) 0183 

X0NE(3)=CKI <ltl)*C0RR(3)+CKI ( 1 , 2 ) *CORR ( 4 ) 0184 

! XONEC 1 )=XONE(3)+XDOTl 0185 

I XONE(2)= (CKK l,l)*CORR(l)+CKI (l,2)»CORR(2) ) 0186 

YONE(2)=CORR(2) 0187 

XINMID=XM0D 0188 

PRINT 559>XINMID ' 0189 

(559 FORMAT! 20H 559 X AT TONE IS ♦F20.10) 0190 

i INTEGRATE FROM T1 TO CAPT 0191 

ELV(1)=BK(1,1) 0192 

i ELV(2)=BK(1,2) 0193 

I ELV(3)=BK(2,1) 0194 

! ELV(4)=BK(2,2) 0195 

' ELV(5)=XMOD ' 0196 

ELV(6)=YMOD 0197 

ELV(7)=0. 0198 

ELV(8)=0. 0199 

DO 302 IA=1,50 • 0200 

DO 303 IC=1,4 0201 

DO 304 ID=1*8 0202 

'304 DE( ID)=ELV( ID)+F( IC)*AK( IC-1 . ID) 0203 

DP( 1 )=-C( 1» 1)*DE( 1 )-C(2»1)*DE(3) 0204 

DP(3)=-C(1,2)*DE(1)-C(2»2)*DE(3) 0205 

DP(2)=-C( 1» I )*DE( 2 )~C( 2»1 )*DE(4) 0206 

DP(4)=-C( 1,2)*DE(2 )-C( 2,2)*DE(4) 0207 

ACON ( 3 ) =DE ( 1 ) *D ( 1 ♦ 1 ) +DE ( 3 ) *D ( 2 » 1 ) + DGUESS* ( DE ( 2 ) *D ( 1 ♦ 1 ) +DE (4)*D(2,1 0208 

+)) 0209 

ACON(4)=DE( 1)*D( 1»2)+DE(3)*D(2»2)+D6UESS*(DE( 2)*D( 1*2)+DE(4)*D( 2,2 0210 

+)) 0211 

SRT=SORTF(ACON( 3)*AC0N(3)+AC0N(4)*AC0N(4) ) 0212 
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CKDET=DE( 1 ) *DE( 4)-DE ( 2 )*DE( 3 ) 

DDET = D( 1»1)*D(2»2)-D( 1 ,2)*D(2,1 ) 

DP( 5)=C( 1 »1 )*DE(5)+C( 1 .2) *DE(6)+(D(1 »1 )*ACON{ 3)+D(l ♦2)*AC0N(4) ) /SR 
+ T 

DP<6)=C(2.1)*OE(5)+C(2»2)*DE(6)+(D(2.1)*ACON(3)+D(2»2)*ACON(4) ) /SR 
+ T 

DP(7)= DDET*DDET*CKDET#CKDET/SRT**3 

DP( 8 )=-DGUESS*DDET*DDET*CK0ET*CKDET/SRT**3 
DO 303 ID=1,8 

303 AK(IC, ID)=TSTEP( 1)*DP( ID) 

IF( IA-1 )801, 801*802 

801 CONTINUE 
XD0T2=DP( 5) 

YD0T2=DP(6) 

802 CONTINUE 

DO 306 ID=1,8 

306 ELV( ID) =ELV( ID) + <Aia 1* I D ) +2 .*AK ( 2 » I D ) +2 .*AK ( 3 » 1 0 ) +AIC ( 4 » I D ) )/6. 

302 T=T+TSTEP(1) 

BKA(1*1)=ELV(1) 

BKA{1*2)=ELV(2) 

BKA(2*1)=ELV(3) 

BKA(2,2)=ELV<4) 

XMOD=ELV( 5) 

YMOD=ELV( 6) 

C0RA=ELV(7) 

C0RB=ELV(8) 

THIS COMPLETES THE INTEGRATION 

XONE { 5 ) =CORR ( 3 ) +BK (1 ♦ 1 ) *( XD0T1-XD0T2 ) +BK ( 2 * 1 ) * ( YDOT 1-YD0T2 ) 

YONEd )=C0RR{4)+BK(1»2 )»(XDOTl-XDOT2)+BK(2*2 )*(YD0T1-YD0T2) 

XBND=XMOD 

YBND=YMOD 

PRINT 558,XBND,YBND 

558 F0RMAT(4H 558»9H X(T)= .F20.10,9H Y(T)» ,F20.10) 

XDIFF=XMID-XINMID 
YDIFF=YFINAL~YBND 

IF( ABSFIXD IFF )+ABSF( YD IFF) -ERROR) 998*998*6001 



0213 

0214'i- 

0215 

0216 
0217 ! 
0218 

0219 

0220 
0221 
0222 

0223 

0224 

0225 

0226 

0227 

0228 

0229 

0230 

0231 

0232 

0233 

0234 

0235 

0236 I 

0237 

0238 

0239 

0240 

0241 

0242 I 

0243 f 

0244 < 

0245 

0246 

0247 

0248 
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li* ■> 



>001 DETER = XONE( 1 ) *CORR ( 2 ) -XONE ( 2 ) *YONE ( 1 ) 

IF ( ABSF ( deter ) - l.E-5 ) 6002 » 6002 *6011 
<•011 DELTAC = XDIFF/X0NE(2) 

DELTAD= ( YDI FF-CORR ( 2 ) *DELTAC ) /CORA 
DELTAT=0. 

>670 IF ( ABSF (DELTAO-. 3) 667 1*6672 *6672 
>672 DELTAC=.5*0ELTAC 
DELTA0=.5*DELTAD 
DELTAT=.5*DELTAT 
60 TO 6670 
>671 GO TO 9001 
>002 DELTAT=XDIFF/X0NE(1) 

DELTAD=(YDIFF-Y0NE(1)*XDIFF/X0NE(1 ) )/CORA 
>673 IF ( ABSF (DELTAD)-. 3 >6674*6675*6675 
>675 DELTAD=.5*DELTAD 
DELTAC=.5*DELTAC 
DELTAT=.5*DELTAT 
GO TO 6673 
>674 CONTINUE 
DELTAC=0. 

DETER = XONE( 1)*C0RA ' 
IF(ABSF(DETER)-1.E-5)9198*9198*9001 

9001 IF ( ABSF (DELTAT >-.5)990 *9002*9002 

9002 DELTAT=.5*DELTAT 
DELTAC=.5*DELTAC 
DELTAD=.5*DELTAD 

I GO TO 9001 

990 PRINT 562*DELTAC*DELTAT*DELTAD 
562 FORMAT! 12H OELTAC = *F20,10*14H DELTA T1 = 

+ = *F20.10> 

TTWO=TONE+DELTAT 
IF(TTW0>996*996*309 
996 T0NE=T0NE*.5 

CGUESS=CGUESS+(TONE/DELTAT*DELTAC> 

DGUESS=DGUESS+(TONE/DELTAT*DElTAD> 

GO TO 310 
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309 IF(TTW0-CAPT)311,995»995 
995 TONE=( TONE+CAPT)*.5 

CGUESS=CGUESS+(CAPT-TONE)*.5/DELTAT*DELTAC 
DGUESS=DGUESS+(CAPT“TONE)*.5/DELTAT*OELTAD 
GO TO 310 
311 TONE=TTWO 

CGUES5=CGUESS+0ELTAC 

DGUESS=DGUESS+OELTAO 

310 CONTINUE 

IF (CAPT-TONE-1. £-6)9771,977 1,9772 

9771 PRINT 9781 

9781 FORMAT! 5H 9781, 16H TONE TOO LARGE) 

GO TO 9991 

9772 IF (TONE-1. £-6)9773,9773, 9774 

9773 PRINT 9782 

9782 FORMAT! 5H 9782, 16H TONE TOO SMALL) 

GO TO 9991 

9774 CONTINUE 
322 CONTINUE 

PRINT 563, CGUESS, TONE 

563 FORMAT !16H 563 NEW C IS ,F20.10,15H NEW TONE IS ,F20.l0) 
999 CONTINUE 
GO TO 1 

998 OTMT=XONE!2)*!CORA*XONE!5)-CORB*YONE!1) ) 

+ +XONE! 1 )*!CORB*CORR !2 )-C0RA*C0RR! 1 ) ) 

IF !ABSF!DTMT)-.01) 9601 ,9601,9602 

9601 IF!ABSF!DTMT)-, 00001)9603, 9603, 9604 

9603 PRINT 9651, XBND 

9651 FORMAT!5H 9651, 9H XMAX* ,F20.10) 

GO TO 9991 

9604 CONTINUE 
GO TO 9602 

9602 EPS2=EPSFAC*DTMT 
DELTAX=XBND-X11 

IF (CELT AX) 8401,840 1,8402 
8401 EPS»EPS*.l 
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GO TO 9627 0321 

J402 EPS=EPS2 0322 

X11=X0ND 0323 

GO TO 9616 0324 

9627 T0NE=T0NE-DELT1T 0325 

DGUESS=DGUESS-DELT1D 0326 

CGUESS=CGUESS-DELT1C 0327 

9616 DELT1C=EPS*(-X0NE( D*CORA) 0328 

DELT1D=EPS*<X0NE( 1 )*CORR( 2 ) -XONE ( 2 ) *YONE ( 1 ) ) 0329 

DELT1T=EPS*(X0NE(2 )*CORA) 0330 

>621 IF(ABSF(DELTlC)+ABSF(DELTlD)+ABSF(DELTlT)-.5)6622.6623»6623 0331 

5623 DELT1C=.5*DELT1C 0332 

DELT1D=.5*DELT1D 0333 

DELT1T=.5*DELT1T 0334 

GO TO 6621 0335 

6622 CONTINUE 0336 

PRINT 9611,DELT1C,DELT1D»DELT1T 0337 

9611 FORMAT(5H 9601»10H DELTAC* ♦F20.10,11H DELTAD= ♦F20.10*11H DEL 0338 
1TAT= ,F20.10) 0339 

CGUESS=CGUESS+DELT1C 0340 

DGUESS=DGUESS+DELT1D 0341 

T0NE=T0NE+DELT1T 0342 

ATANC=ATANF(CGUESS ) 0343 

ATAND=ATANF(DGUESS) 0344 

PRINT 9617,ATANC»ATAND 0345 

9617 FORMAT(5H 9617*6H C= ,F20.10*10H RADIANS* *6H D= ,F20.10»9H R 0346 

; +ADIANS) 0347 

DELTC=ATANC-ATC11 0348 

DELTD=ATAND-AT011 0349 

IF( ABSF(DELTC)+ABSF(DELTD)-l.E-6) 9981 *9981 »9982 0350 

9981 PRINT 9983 0351 

9983 FORMAT (5H 9983*49H CORRECTION IS LESS THAN ,000001 - WE ARE DON 0352 

IE//) 0353 

GO TO 9991 ' 0354 

9982^ATC11=ATANC \ 0355 

ATD11=ATAND 0356 
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9985 CONTINUE 

PRINT 9612»EPS ^350 

9612 FORMAT! 12H 9612 EPS= »F20.10) 0359 

AAK(1*1)=BK{1.1)+CGUESS#BK(1»2) 0360 

AAK(2*1)=BK(2»1)+CGUESS*BK(2»2) 0361 

AAK ( 1 ♦ 2 ) =BK (1.1) +DGUESS*BK (1.2) 0362 

AAK(2.2)=BK(2.1)+OGUESS*BK(2.2) 0363 

PRINT 9987 036A 

9987 FORMAT(5H 9987.18H K* AT TONE M I NUS . lOX ♦ 1 5HK* AT TONE PLUS//) 0365 

PRINT 9988. ( (AAK( I »J) »J*1»2) .1=1.2) 0366 

9988 FORMAT! 10X.F12.5.14X.F12.5) 0367 

PRINT 9989 0360 

9989 FORMAT! IHO. 12H XOOT MI NUS . lOX . 9HXDOT PLUS) 0369 

PRINT 9990.XDOT1.XDOT2.YDOT1.YDOT2 0370 

9990 FORMAT(5X.F12.5.10X.F12.5) 0371 

GO TO 1 0372 

9198 PRINT 9199 0373 

9199 FORMAT!68H RAN)C OF FIRST CORRECTION MATRIX IS ONE - NO ADMISSIBLE 0374 

. ICURVES EXIST) 0375 

9991 CONTINUE 0376 

STOP 0377 

END 0370 

END 0379 
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